Info from Kim:
library(igraph)
##
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
## The following object is masked from 'package:base':
##
## union
library(rgl)
library(rafalib)
library(slingshot)
## Loading required package: princurve
library(msigdbr)
library(fgsea)
## Loading required package: Rcpp
library(knitr)
library(leiden)
sds <- readRDS("../data/finalTrajectory/sling.rds")
clDatta <- readRDS("../data/finalTrajectory/dattaCl_noResp_noMV.rds")
cols <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#FFFF33", "#A65628", "#F781BF", "#999999")
view3d( theta = 10, phi = 10)
# Datta labels
plot3d(reducedDim(sds), aspect = 'iso', col=cols[clDatta], alpha=.6, xlab="UMAP1", ylab="UMAP2", zlab="UMAP3", box=FALSE, top=TRUE)
load("../data/list_TF_Kim.RData")
HBC* - GBC - iOSN - mOSN
## activated HBC
load("../data/associationNetworkKim/result_realdata/3-dataTFtranscelsalpha059515result.RData")
adjHBCAct <- adj.zinb1
dimnames(adjHBCAct) <- list(colnames(datatest), colnames(datatest))
adjHBCAct <- adjHBCAct[!rowSums(adjHBCAct) == 0,]
adjHBCAct <- adjHBCAct[,!colSums(adjHBCAct) == 0]
rm(adj.zinb1, datatest)
## GBC
load("../data/associationNetworkKim/result_realdata/4-dataTFtranscelsalpha059515result.RData")
adjGBC <- adj.zinb1
dimnames(adjGBC) <- list(colnames(datatest), colnames(datatest))
adjGBC <- adjGBC[!rowSums(adjGBC) == 0,]
adjGBC <- adjGBC[,!colSums(adjGBC) == 0]
rm(adj.zinb1, datatest)
## iOSN
load("../data/associationNetworkKim/result_realdata/2-dataTFtranscelsalpha059515result.RData")
adjiOSN <- adj.zinb1
dimnames(adjiOSN) <- list(colnames(datatest), colnames(datatest))
adjiOSN <- adjiOSN[!rowSums(adjiOSN) == 0,]
adjiOSN <- adjiOSN[,!colSums(adjiOSN) == 0]
rm(adj.zinb1, datatest)
## mOSN
load("../data/associationNetworkKim/result_realdata/5-dataTFtranscelsalpha059515result.RData")
adjmOSN <- adj.zinb1
dimnames(adjmOSN) <- list(colnames(datatest), colnames(datatest))
adjmOSN <- adjmOSN[!rowSums(adjmOSN) == 0,]
adjmOSN <- adjmOSN[,!colSums(adjmOSN) == 0]
rm(adj.zinb1, datatest)
# gHBC <- graph_from_adjacency_matrix(adjHBC,mode="undirected")
# cols <- c(brewer.pal(8, "Dark2"),brewer.pal(8, "Set2"))
# V(g)$color <- cols[cl$cluster[order(cl$X, decreasing=FALSE)]+2]
# lo <- layout_with_fr(gHBC)
# plot(g,vertex.size=3,vertex.label=NA,vertex.frame.color=V(g)$color,layout=lo,edge.width=0.6)
# legend("bottomleft",paste("comm",seq(1,6)),col=cols[3:8],pch=16, bty='n')
degreeDf <- data.frame(degree=1:52,
hbcAct=NA,
gbc=NA,
iOSN=NA,
mOSN=NA)
degreeHBCAct <- table(rowSums(adjHBCAct))
degreeDf$hbcAct[match(names(degreeHBCAct), degreeDf$degree)] <- degreeHBCAct
degreeGBC <- table(rowSums(adjGBC))
degreeDf$gbc[match(names(degreeGBC), degreeDf$degree)] <- degreeGBC
degreeiOSN <- table(rowSums(adjiOSN))
degreeDf$iOSN[match(names(degreeiOSN), degreeDf$degree)] <- degreeiOSN
degreemOSN <- table(rowSums(adjmOSN))
degreeDf$mOSN[match(names(degreemOSN), degreeDf$degree)] <- degreemOSN
mypar(mfrow=c(2,2))
barplot(degreeDf$hbcAct, names=degreeDf$degree, main="HBC*: Node-level degree", ylim=c(0,400))
barplot(degreeDf$gbc, names=degreeDf$degree, main="GBC: Node-level degree", ylim=c(0,400))
barplot(degreeDf$iOSN, names=degreeDf$degree, main="iOSN: Node-level degree", ylim=c(0,400))
barplot(degreeDf$mOSN, names=degreeDf$degree, main="mOSN: Node-level degree", ylim=c(0,400))
extractCore <- function(A, degree = 3){
converg <- FALSE
old.nrow <- nrow(A)
while(!converg){
d <- colSums(A)
to.keep <- which(d>=degree)
if(old.nrow==length(to.keep)){
converg <- TRUE
}
old.nrow <- length(to.keep)
A <- A[to.keep,to.keep]
}
return(A)
}
adjHBCAct2 <- extractCore(adjHBCAct, degree = 2)
adjGBC2 <- extractCore(adjGBC, degree = 2)
adjiOSN2 <- extractCore(adjiOSN, degree = 2)
adjmOSN2 <- extractCore(adjmOSN, degree = 2)
mypar(mfrow=c(2,2))
gHBCAct2 <- graph_from_adjacency_matrix(adjHBCAct2,mode="undirected")
plot(gHBCAct2,vertex.size=3,vertex.label=NA,edge.width=0.6, main="HBC*")
gGBC2 <- graph_from_adjacency_matrix(adjGBC2,mode="undirected")
plot(gGBC2,vertex.size=3,vertex.label=NA,edge.width=0.6, main="GBC")
giOSN2 <- graph_from_adjacency_matrix(adjiOSN2,mode="undirected")
plot(gGBC2,vertex.size=3,vertex.label=NA,edge.width=0.6, main="iOSN")
gmOSN2 <- graph_from_adjacency_matrix(adjmOSN2,mode="undirected")
plot(gGBC2,vertex.size=3,vertex.label=NA,edge.width=0.6, main="mOSN")
saveRDS(list("hbcAct" = adjHBCAct2,
"GBC" = adjGBC2,
"iOSN" = adjiOSN2,
"mOSN" = adjmOSN2), file="../data/adjMatrices_2core.rds")
library(RColorBrewer)
library(HCD)
gradualPlot <- function(cl, A){
ncl <- length(unique(cl))
for(cc in 1:ncl){
id <- which(cl %in% (1:cc))
curG <- graph_from_adjacency_matrix(A[id,id],mode="undirected")
V(curG)$color <- cols[cl[names(V(curG))]]
lo <- layout_with_fr(curG)
plot(curG,
vertex.size=3,
vertex.label=NA,
edge.width=0.6,
vertex.frame.color=cols[cl[id]],
color=V(curG)$color,
layout = lo
)
}
}
## C5 category is according to gene ontology grouping: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4707969/pdf/nihms-743907.pdf
geneSets <- msigdbr(species = "Mus musculus", category = "C5", subcategory = "BP")
gsea <- function(genes, background, geneSets, n=10, minSize=5, name=NULL){
### filter background to only include genes that we assessed.
geneSets <- geneSets[geneSets$gene_symbol %in% background,]
m_list <- geneSets %>% split(x = .$gene_symbol, f = .$gs_name)
# gene set must have at least minSize genes in background.
m_list <- m_list[unlist(lapply(m_list, length)) >= minSize]
overlapPval <- unlist(lapply(m_list, function(gs){
# genes in community and gene set
inBoth <- sum(genes %in% gs)
# genes in community and not in gene set
inComOnly <- length(genes) - inBoth
# genes in background and gene set
inGsBack <- sum(background %in% gs)
# genes in background and not in gene set
outGsBack <- length(background) - inGsBack
m <- matrix(c(inBoth, inComOnly,
inGsBack, outGsBack),
nrow =2, ncol=2, byrow=TRUE,
dimnames = list(c("in community", "out community"),
c("in gene set", "out gene set")))
fis <- fisher.test(m, alternative = "greater")
pval <- fis$p.value
return(pval)
}))
padj <- p.adjust(overlapPval, "fdr")
oo <- order(overlapPval, decreasing=FALSE)
res <- data.frame(geneSet = names(m_list)[oo[1:n]],
pval = overlapPval[oo[1:n]],
padj = padj[oo[1:n]],
row.names = NULL)
kable(res, caption=name, label=name)
}
#set.seed(2)
clHBCAct <- leiden(adjHBCAct2, resolution_parameter=.4, seed = 7)
table(clHBCAct)
## clHBCAct
## 1 2 3 4
## 378 286 234 28
gHBCAct <- graph_from_adjacency_matrix(adjHBCAct2,mode="undirected")
cols <- brewer.pal(8, "Dark2")
V(gHBCAct)$color <- cols[clHBCAct[names(V(gHBCAct))]]
plot(gHBCAct,
vertex.size=3,
vertex.label=NA,
#vertex.frame.color=V(gHBC)$color,
color=V(gHBCAct)$color,
edge.width=0.6)
comHBCAct <- make_clusters(gHBCAct, membership = clHBCAct)
plot(comHBCAct, gHBCAct,vertex.size=3, vertex.label=NA, edge.width=0.6)
## gradual plot
gradualPlot(clHBCAct, adjHBCAct2)
## gene set enrichment on HCD results
ncl <- length(unique(comHBCAct$membership))
kab <- list()
# kabNames <- c("Krt / cell cycle / stimulus",
# "Chromatin",
# "Epigenetics, cell cycle",
# "Cell cycle")
kabNames <- c("Cell cycle, DNA damage",
"Epigenetics",
"(Epithelial) cell differentiation",
"Cell cycle, DNA replication",
"DNA replication")
for(kk in 1:ncl){
genes <- rownames(adjHBCAct2)[which(comHBCAct$membership == kk)]
kab[[kk]] <- gsea(genes = genes,
background = TF.list,
geneSets = geneSets,
name = kabNames[kk])
}
kab
## [[1]]
##
##
## Table: Cell cycle, DNA damage
##
## geneSet pval padj
## --------------------------------------------- ---------- ----------
## GO_CELLULAR_RESPONSE_TO_DNA_DAMAGE_STIMULUS 0.0010309 0.6004173
## GO_REGULATION_OF_CELL_CYCLE 0.0023338 0.6004173
## GO_CARBOHYDRATE_HOMEOSTASIS 0.0023584 0.6004173
## GO_RESPONSE_TO_CARBOHYDRATE 0.0029978 0.6004173
## GO_PROCESS_UTILIZING_AUTOPHAGIC_MECHANISM 0.0033514 0.6004173
## GO_REGULATION_OF_MITOTIC_CELL_CYCLE 0.0034382 0.6004173
## GO_NEGATIVE_REGULATION_OF_CELL_CYCLE_PROCESS 0.0037266 0.6004173
## GO_REGULATION_OF_AUTOPHAGY 0.0041650 0.6004173
## GO_DNA_REPAIR 0.0044003 0.6004173
## GO_REGULATION_OF_CELL_CYCLE_PHASE_TRANSITION 0.0045432 0.6004173
##
## [[2]]
##
##
## Table: Epigenetics
##
## geneSet pval padj
## -------------------------------------------------------- ---------- -----
## GO_PROTEIN_DEMETHYLATION 0.0137397 1
## GO_DEMETHYLATION 0.0196971 1
## GO_REGULATION_OF_DNA_REPLICATION 0.0275424 1
## GO_FOAM_CELL_DIFFERENTIATION 0.0278109 1
## GO_HISTONE_H3_K4_DEMETHYLATION 0.0385626 1
## GO_HISTONE_H3_K4_DEMETHYLATION_TRIMETHYL_H3_K4_SPECIFIC 0.0385626 1
## GO_HISTONE_H3_K4_METHYLATION 0.0479833 1
## GO_RIBONUCLEOPROTEIN_COMPLEX_SUBUNIT_ORGANIZATION 0.0550234 1
## GO_REGULATION_OF_CHOLESTEROL_BIOSYNTHETIC_PROCESS 0.0708770 1
## GO_STEROL_BIOSYNTHETIC_PROCESS 0.0708770 1
##
## [[3]]
##
##
## Table: (Epithelial) cell differentiation
##
## geneSet pval padj
## ------------------------------------------------ ---------- ----------
## GO_SCHWANN_CELL_DIFFERENTIATION 0.0037526 0.9447078
## GO_REGULATION_OF_CELL_POPULATION_PROLIFERATION 0.0051764 0.9447078
## GO_EPIDERMIS_DEVELOPMENT 0.0064884 0.9447078
## GO_POSITIVE_REGULATION_OF_CELL_DEATH 0.0072210 0.9447078
## GO_PERIPHERAL_NERVOUS_SYSTEM_DEVELOPMENT 0.0106454 0.9447078
## GO_RAS_PROTEIN_SIGNAL_TRANSDUCTION 0.0117230 0.9447078
## GO_EPITHELIUM_DEVELOPMENT 0.0134552 0.9447078
## GO_NEGATIVE_REGULATION_OF_BMP_SIGNALING_PATHWAY 0.0137142 0.9447078
## GO_REGULATION_OF_CELL_CYCLE 0.0159675 0.9447078
## GO_APOPTOTIC_PROCESS 0.0161161 0.9447078
##
## [[4]]
##
##
## Table: Cell cycle, DNA replication
##
## geneSet pval padj
## --------------------------------------------- ---------- ----------
## GO_DNA_REPLICATION_INITIATION 0.0000000 0.0001121
## GO_DNA_DEPENDENT_DNA_REPLICATION 0.0000081 0.0075936
## GO_CELL_CYCLE_G1_S_PHASE_TRANSITION 0.0000101 0.0075936
## GO_DNA_UNWINDING_INVOLVED_IN_DNA_REPLICATION 0.0000158 0.0089155
## GO_DNA_GEOMETRIC_CHANGE 0.0001095 0.0494663
## GO_CELL_CYCLE_PHASE_TRANSITION 0.0002037 0.0766622
## GO_DNA_REPLICATION 0.0003546 0.1041845
## GO_MITOTIC_CELL_CYCLE 0.0003691 0.1041845
## GO_DNA_METABOLIC_PROCESS 0.0009372 0.2351229
## GO_DNA_CONFORMATION_CHANGE 0.0013729 0.3100113
set.seed(3)
clGBC <- leiden(adjGBC2, resolution_parameter=.4, seed=7)
table(clGBC)
## clGBC
## 1 2 3 4
## 331 279 274 107
gGBC <- graph_from_adjacency_matrix(adjGBC2,mode="undirected")
cols <- brewer.pal(8, "Dark2")
V(gGBC)$color <- cols[clGBC[names(V(gGBC))]]
plot(gGBC,
vertex.size=3,
vertex.label=NA,
#vertex.frame.color=V(gHBC)$color,
color=V(gGBC)$color,
edge.width=0.6)
comGBC <- make_clusters(gGBC, membership = clGBC)
plot(comGBC, gGBC,vertex.size=3, vertex.label=NA, edge.width=0.6)
## gradual plot
gradualPlot(clGBC, adjGBC2)
## gene set enrichment
ncl <- length(unique(comGBC$membership))
kab <- list()
# kabNames <- c("Chromatin, neuron/cell development",
# "P53, DNA replication, cell cycle",
# "Cell differentiation / signaling",
# "Notch, signaling",
# "DNA repair")
kabNames <- c("P53, DNA replication",
"Notch signaling",
"Cell differentiation / signaling",
"RNA splicing, DNA repair, epigenetics",
"NIK/NF Kappa-B signaling")
for(kk in 1:ncl){
genes <- rownames(adjGBC2)[which(comGBC$membership == kk)]
print(gsea(genes = genes,
background = TF.list,
geneSets = geneSets,
name = kabNames[kk]))
}
##
##
## Table: P53, DNA replication
##
## geneSet pval padj
## ------------------------------------------------------------ ---------- -----
## GO_REGULATION_OF_SIGNAL_TRANSDUCTION_BY_P53_CLASS_MEDIATOR 0.0073943 1
## GO_DNA_REPLICATION_INITIATION 0.0088812 1
## GO_HISTONE_PHOSPHORYLATION 0.0088812 1
## GO_SIGNAL_TRANSDUCTION_BY_P53_CLASS_MEDIATOR 0.0111609 1
## GO_DNA_REPLICATION 0.0173178 1
## GO_SNRNA_TRANSCRIPTION 0.0176173 1
## GO_NCRNA_TRANSCRIPTION 0.0184906 1
## GO_TRANSCRIPTION_ELONGATION_FROM_RNA_POLYMERASE_II_PROMOTER 0.0212404 1
## GO_POSITIVE_REGULATION_OF_CHROMOSOME_ORGANIZATION 0.0277202 1
## GO_DNA_TEMPLATED_TRANSCRIPTION_ELONGATION 0.0301290 1
##
##
## Table: Notch signaling
##
## geneSet pval padj
## ----------------------------------------------- ---------- ----------
## GO_REGULATION_OF_CELL_POPULATION_PROLIFERATION 0.0003936 0.3441872
## GO_REGULATION_OF_NOTCH_SIGNALING_PATHWAY 0.0004437 0.3441872
## GO_RESPONSE_TO_NITROGEN_COMPOUND 0.0005204 0.3441872
## GO_NEGATIVE_REGULATION_OF_SIGNALING 0.0007461 0.3441872
## GO_CELLULAR_RESPONSE_TO_EXTERNAL_STIMULUS 0.0010464 0.3441872
## GO_RESPONSE_TO_PEPTIDE 0.0013035 0.3441872
## GO_RESPONSE_TO_STARVATION 0.0016670 0.3441872
## GO_CELLULAR_RESPONSE_TO_STARVATION 0.0016766 0.3441872
## GO_CARDIOVASCULAR_SYSTEM_DEVELOPMENT 0.0017096 0.3441872
## GO_NOTCH_SIGNALING_PATHWAY 0.0020545 0.3441872
##
##
## Table: Cell differentiation / signaling
##
## geneSet pval padj
## -------------------------------------------------- ---------- -----
## GO_RNA_SPLICING_VIA_TRANSESTERIFICATION_REACTIONS 0.0055550 1
## GO_REGULATION_OF_GENE_SILENCING 0.0070297 1
## GO_RNA_SPLICING 0.0077460 1
## GO_CELLULAR_PROTEIN_CONTAINING_COMPLEX_ASSEMBLY 0.0084200 1
## GO_DNA_REPAIR 0.0088367 1
## GO_MRNA_METABOLIC_PROCESS 0.0104386 1
## GO_TRANSCRIPTION_PREINITIATION_COMPLEX_ASSEMBLY 0.0106774 1
## GO_REGULATION_OF_GENE_EXPRESSION_EPIGENETIC 0.0139632 1
## GO_REGULATION_OF_STEM_CELL_DIFFERENTIATION 0.0142288 1
## GO_REGULATION_OF_DNA_BINDING 0.0146403 1
##
##
## Table: RNA splicing, DNA repair, epigenetics
##
## geneSet pval padj
## ------------------------------------------------------------- ---------- -----
## GO_REGULATION_OF_NIK_NF_KAPPAB_SIGNALING 0.0034495 1
## GO_ORGANELLE_LOCALIZATION 0.0051168 1
## GO_NIK_NF_KAPPAB_SIGNALING 0.0164556 1
## GO_ESTABLISHMENT_OF_ORGANELLE_LOCALIZATION 0.0206255 1
## GO_MICROTUBULE_CYTOSKELETON_ORGANIZATION_INVOLVED_IN_MITOSIS 0.0227626 1
## GO_POSITIVE_REGULATION_OF_NIK_NF_KAPPAB_SIGNALING 0.0298471 1
## GO_SPINDLE_ASSEMBLY 0.0298471 1
## GO_HISTONE_H3_ACETYLATION 0.0309471 1
## GO_HISTONE_H4_ACETYLATION 0.0355192 1
## GO_CELLULAR_RESPONSE_TO_CALCIUM_ION 0.0379565 1
kab
## list()
set.seed(4)
cliOSN <- leiden(adjiOSN2, resolution_parameter=.3, seed=7)
table(cliOSN)
## cliOSN
## 1 2 3 4 5 6 7
## 235 83 77 27 21 16 6
giOSN <- graph_from_adjacency_matrix(adjiOSN2,mode="undirected")
cols <- brewer.pal(8, "Dark2")
V(giOSN)$color <- cols[cliOSN[names(V(giOSN))]]
plot(giOSN,
vertex.size=3,
vertex.label=NA,
color=V(giOSN)$color,
edge.width=0.6)
# igraph plot
comiOSN <- make_clusters(giOSN, membership = cliOSN)
plot(comiOSN, giOSN,vertex.size=3, vertex.label=NA, edge.width=0.6)
## gradual plot
gradualPlot(cliOSN, adjiOSN2)
## gene set enrichment on HCD results
ncl <- length(unique(comiOSN$membership))
kab <- list()
kabNames <- c("Neurogenesis",
"Defense (virus) response",
"P53, membrane",
"Axon extension, cytokinesis, epigenetic gene expression regulation",
"Cytokinesis, DNA replication",
"Interleukins, cell cycle, protein localization",
"Stem cell population maintenance, transcription")
for(kk in 1:ncl){
genes <- rownames(adjiOSN2)[which(comiOSN$membership == kk)]
print(gsea(genes = genes,
background = TF.list,
geneSets = geneSets,
name = kabNames[kk]))
}
##
##
## Table: Neurogenesis
##
## geneSet pval padj
## -------------------------------------------------------- ---------- ----------
## GO_TELENCEPHALON_DEVELOPMENT 0.0002684 0.3249968
## GO_NEUROGENESIS 0.0002879 0.3249968
## GO_NEURON_DIFFERENTIATION 0.0008345 0.3692460
## GO_REGULATION_OF_NERVOUS_SYSTEM_DEVELOPMENT 0.0008941 0.3692460
## GO_EMBRYO_DEVELOPMENT 0.0008958 0.3692460
## GO_CENTRAL_NERVOUS_SYSTEM_DEVELOPMENT 0.0011003 0.3692460
## GO_POSITIVE_REGULATION_OF_CELL_POPULATION_PROLIFERATION 0.0011447 0.3692460
## GO_FOREBRAIN_DEVELOPMENT 0.0014588 0.3837962
## GO_NEGATIVE_REGULATION_OF_NEURON_DIFFERENTIATION 0.0015983 0.3837962
## GO_REGIONALIZATION 0.0019554 0.3837962
##
##
## Table: Defense (virus) response
##
## geneSet pval padj
## ------------------------------------------------------------------- ---------- -----
## GO_DEFENSE_RESPONSE_TO_VIRUS 0.0032635 1
## GO_RESPONSE_TO_VIRUS 0.0038295 1
## GO_RESPONSE_TO_BIOTIC_STIMULUS 0.0056459 1
## GO_NEGATIVE_REGULATION_OF_JNK_CASCADE 0.0118638 1
## GO_REGULATION_OF_MRNA_CATABOLIC_PROCESS 0.0119217 1
## GO_REGULATION_OF_STRESS_ACTIVATED_PROTEIN_KINASE_SIGNALING_CASCADE 0.0119217 1
## GO_RESPONSE_TO_TYPE_I_INTERFERON 0.0135956 1
## GO_CELLULAR_HOMEOSTASIS 0.0159622 1
## GO_DEFENSE_RESPONSE_TO_OTHER_ORGANISM 0.0176323 1
## GO_RESPONSE_TO_CYTOKINE 0.0185354 1
##
##
## Table: P53, membrane
##
## geneSet pval padj
## ----------------------------------------------------------------------------------------- ---------- -----
## GO_MITOCHONDRIAL_MEMBRANE_ORGANIZATION 0.0209248 1
## GO_PROTEIN_MONOUBIQUITINATION 0.0209248 1
## GO_REGULATION_OF_MEMBRANE_PERMEABILITY 0.0209248 1
## GO_INTRINSIC_APOPTOTIC_SIGNALING_PATHWAY_BY_P53_CLASS_MEDIATOR 0.0242866 1
## GO_INTRINSIC_APOPTOTIC_SIGNALING_PATHWAY_IN_RESPONSE_TO_DNA_DAMAGE_BY_P53_CLASS_MEDIATOR 0.0310897 1
## GO_MAINTENANCE_OF_PROTEIN_LOCATION 0.0310897 1
## GO_REGULATION_OF_GENE_SILENCING 0.0312073 1
## GO_MRNA_TRANSCRIPTION 0.0369826 1
## GO_REGULATION_OF_CHROMATIN_SILENCING 0.0369826 1
## GO_BASE_EXCISION_REPAIR 0.0400952 1
##
##
## Table: Axon extension, cytokinesis, epigenetic gene expression regulation
##
## geneSet pval padj
## ------------------------------------------------- ---------- -----
## GO_REGULATION_OF_VESICLE_MEDIATED_TRANSPORT 0.0718594 1
## GO_DNA_TEMPLATED_TRANSCRIPTION_TERMINATION 0.0870021 1
## GO_REGULATION_OF_GENE_EXPRESSION_EPIGENETIC 0.0879146 1
## GO_ATF6_MEDIATED_UNFOLDED_PROTEIN_RESPONSE 0.0989998 1
## GO_AXON_EXTENSION 0.0989998 1
## GO_CHOLESTEROL_EFFLUX 0.0989998 1
## GO_CHOLESTEROL_STORAGE 0.0989998 1
## GO_CYTOKINESIS 0.0989998 1
## GO_MITOTIC_SPINDLE_ASSEMBLY 0.0989998 1
## GO_NEGATIVE_REGULATION_OF_INNATE_IMMUNE_RESPONSE 0.0989998 1
##
##
## Table: Cytokinesis, DNA replication
##
## geneSet pval padj
## ------------------------------------------------------------------- ---------- -----
## GO_CYTOKINESIS 0.0034643 1
## GO_POSITIVE_REGULATION_OF_DNA_DEPENDENT_DNA_REPLICATION 0.0034643 1
## GO_REGULATION_OF_CYTOKINESIS 0.0034643 1
## GO_TROPHOBLAST_GIANT_CELL_DIFFERENTIATION 0.0034643 1
## GO_EMBRYO_DEVELOPMENT_ENDING_IN_BIRTH_OR_EGG_HATCHING 0.0062909 1
## GO_CELL_DIFFERENTIATION_INVOLVED_IN_EMBRYONIC_PLACENTA_DEVELOPMENT 0.0072449 1
## GO_POSITIVE_REGULATION_OF_DNA_REPLICATION 0.0072449 1
## GO_POSITIVE_REGULATION_OF_CELL_CYCLE_ARREST 0.0087899 1
## GO_EMBRYONIC_PLACENTA_DEVELOPMENT 0.0095574 1
## GO_CELL_CYCLE_DNA_REPLICATION 0.0104550 1
##
##
## Table: Interleukins, cell cycle, protein localization
##
## geneSet pval padj
## ---------------------------------------------------------------- ---------- -----
## GO_RESPONSE_TO_INTERLEUKIN_4 0.0072146 1
## GO_G0_TO_G1_TRANSITION 0.0168525 1
## GO_CIRCADIAN_RHYTHM 0.0284399 1
## GO_POSITIVE_REGULATION_OF_CELLULAR_PROTEIN_LOCALIZATION 0.0341263 1
## GO_DNA_INTEGRITY_CHECKPOINT 0.0482564 1
## GO_RHYTHMIC_PROCESS 0.0508153 1
## GO_POSITIVE_REGULATION_OF_ESTABLISHMENT_OF_PROTEIN_LOCALIZATION 0.0560070 1
## GO_ESTABLISHMENT_OF_PROTEIN_LOCALIZATION_TO_ORGANELLE 0.0586845 1
## GO_CHROMATIN_SILENCING_AT_TELOMERE 0.0601135 1
## GO_HISTONE_H3_K9_ACETYLATION 0.0601135 1
##
##
## Table: Stem cell population maintenance, transcription
##
## geneSet pval padj
## -------------------------------------------------------------------- ---------- -----
## GO_REGULATION_OF_STEM_CELL_POPULATION_MAINTENANCE 0.0306392 1
## GO_REGULATION_OF_DNA_TEMPLATED_TRANSCRIPTION_ELONGATION 0.0456629 1
## GO_REGULATION_OF_PROTEIN_COMPLEX_ASSEMBLY 0.1352331 1
## GO_REGULATION_OF_CELLULAR_AMIDE_METABOLIC_PROCESS 0.1386647 1
## GO_PEPTIDE_BIOSYNTHETIC_PROCESS 0.1522779 1
## GO_AMIDE_BIOSYNTHETIC_PROCESS 0.1623697 1
## GO_DNA_TEMPLATED_TRANSCRIPTION_ELONGATION 0.1623697 1
## GO_POSITIVE_REGULATION_OF_DNA_BINDING_TRANSCRIPTION_FACTOR_ACTIVITY 0.1623697 1
## GO_CELLULAR_AMIDE_METABOLIC_PROCESS 0.1756696 1
## GO_MAINTENANCE_OF_CELL_NUMBER 0.1952891 1
kab
## list()
set.seed(2)
clmOSN <- leiden(adjmOSN2, resolution_parameter=.4, seed=7)
table(clmOSN)
## clmOSN
## 1 2 3 4 5
## 329 182 139 81 41
gmOSN <- graph_from_adjacency_matrix(adjmOSN2,mode="undirected")
cols <- brewer.pal(8, "Dark2")
V(gmOSN)$color <- cols[clmOSN[names(V(gmOSN))]]
plot(gmOSN,
vertex.size=3,
vertex.label=NA,
color=V(gmOSN)$color,
edge.width=0.6)
# igraph plot
commOSN <- make_clusters(gmOSN, membership = clmOSN)
plot(commOSN, gmOSN,vertex.size=3, vertex.label=NA, edge.width=0.6)
## gradual plot
gradualPlot(clmOSN, adjmOSN2)
## gene set enrichment on HCD results
ncl <- length(unique(commOSN$membership))
kab <- list()
kabNames <- c("Chromatin",
"Chromatin, transcription",
"TGBF-Beta pathway, neuron death",
"Viral processes",
"Cell projection, others")
for(kk in 1:ncl){
genes <- rownames(adjmOSN2)[which(commOSN$membership == kk)]
print(gsea(genes = genes,
background = TF.list,
geneSets = geneSets,
name = kabNames[kk]))
}
##
##
## Table: Chromatin
##
## geneSet pval padj
## -------------------------------------------------- ---------- -----
## GO_COVALENT_CHROMATIN_MODIFICATION 0.0076200 1
## GO_CHROMATIN_ORGANIZATION 0.0088562 1
## GO_PROTEIN_ACYLATION 0.0090561 1
## GO_CHROMOSOME_ORGANIZATION 0.0124554 1
## GO_RNA_SPLICING_VIA_TRANSESTERIFICATION_REACTIONS 0.0151444 1
## GO_RRNA_TRANSCRIPTION 0.0160161 1
## GO_POSITIVE_REGULATION_OF_DEFENSE_RESPONSE 0.0176331 1
## GO_PROTEIN_ACETYLATION 0.0179776 1
## GO_PEPTIDYL_LYSINE_ACETYLATION 0.0196668 1
## GO_REGULATION_OF_RNA_SPLICING 0.0210600 1
##
##
## Table: Chromatin, transcription
##
## geneSet pval padj
## ----------------------------------------------------------------------- ---------- -----
## GO_ATP_DEPENDENT_CHROMATIN_REMODELING 0.0099950 1
## GO_DNA_TEMPLATED_TRANSCRIPTION_ELONGATION 0.0113455 1
## GO_PROTEIN_DNA_COMPLEX_SUBUNIT_ORGANIZATION 0.0202977 1
## GO_FERTILIZATION 0.0218270 1
## GO_RESPONSE_TO_EXTRACELLULAR_STIMULUS 0.0219166 1
## GO_RESPONSE_TO_OSMOTIC_STRESS 0.0300984 1
## GO_CHROMATIN_ORGANIZATION 0.0313435 1
## GO_INSULIN_RECEPTOR_SIGNALING_PATHWAY 0.0338212 1
## GO_ORGANOPHOSPHATE_BIOSYNTHETIC_PROCESS 0.0338212 1
## GO_CELLULAR_PROCESS_INVOLVED_IN_REPRODUCTION_IN_MULTICELLULAR_ORGANISM 0.0358365 1
##
##
## Table: TGBF-Beta pathway, neuron death
##
## geneSet pval padj
## ---------------------------------------------------------------------------- ---------- ----------
## GO_TRANSFORMING_GROWTH_FACTOR_BETA_RECEPTOR_SIGNALING_PATHWAY 0.0001683 0.3800204
## GO_CELLULAR_RESPONSE_TO_INORGANIC_SUBSTANCE 0.0005379 0.6072730
## GO_RESPONSE_TO_TRANSFORMING_GROWTH_FACTOR_BETA 0.0020270 0.8644903
## GO_REGULATION_OF_DNA_BINDING 0.0021189 0.8644903
## GO_NEURON_DEATH 0.0027505 0.8644903
## GO_REGULATION_OF_DNA_TEMPLATED_TRANSCRIPTION_IN_RESPONSE_TO_STRESS 0.0027646 0.8644903
## GO_TRANSMEMBRANE_RECEPTOR_PROTEIN_SERINE_THREONINE_KINASE_SIGNALING_PATHWAY 0.0028099 0.8644903
## GO_POSITIVE_REGULATION_OF_NEURON_DEATH 0.0030629 0.8644903
## GO_CELLULAR_RESPONSE_TO_CALCIUM_ION 0.0038982 0.9780268
## GO_RESPONSE_TO_METAL_ION 0.0049149 1.0000000
##
##
## Table: Viral processes
##
## geneSet pval padj
## ------------------------------------------------------------------------------------------ ---------- -----
## GO_VIRAL_LIFE_CYCLE 0.0010567 1
## GO_INTERSPECIES_INTERACTION_BETWEEN_ORGANISMS 0.0016139 1
## GO_PROTEIN_MONOUBIQUITINATION 0.0039224 1
## GO_REGULATION_OF_VIRAL_LIFE_CYCLE 0.0051478 1
## GO_PROTEIN_MODIFICATION_BY_SMALL_PROTEIN_CONJUGATION 0.0060424 1
## GO_REGULATION_OF_NUCLEAR_TRANSCRIBED_MRNA_CATABOLIC_PROCESS_DEADENYLATION_DEPENDENT_DECAY 0.0080769 1
## GO_REGULATION_OF_SYMBIOSIS_ENCOMPASSING_MUTUALISM_THROUGH_PARASITISM 0.0103493 1
## GO_PROTEIN_SUMOYLATION 0.0125485 1
## GO_POSITIVE_REGULATION_OF_VIRAL_PROCESS 0.0166502 1
## GO_CELLULAR_MACROMOLECULE_CATABOLIC_PROCESS 0.0180804 1
##
##
## Table: Cell projection, others
##
## geneSet pval padj
## ----------------------------------------------------------- ---------- -----
## GO_CELL_PROJECTION_ASSEMBLY 0.0135173 1
## GO_REGULATION_OF_CIRCADIAN_RHYTHM 0.0137572 1
## GO_NEGATIVE_REGULATION_OF_SMOOTHENED_SIGNALING_PATHWAY 0.0165868 1
## GO_ENDOSOMAL_TRANSPORT 0.0209788 1
## GO_POSITIVE_REGULATION_OF_CELL_CYCLE_G2_M_PHASE_TRANSITION 0.0209788 1
## GO_SPINAL_CORD_PATTERNING 0.0209788 1
## GO_PROTEIN_DEPHOSPHORYLATION 0.0257976 1
## GO_PROTEIN_LOCALIZATION_TO_MEMBRANE 0.0437886 1
## GO_CELL_CYCLE_G2_M_PHASE_TRANSITION 0.0474675 1
## GO_MITOTIC_SISTER_CHROMATID_SEGREGATION 0.0488742 1
kab
## list()